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Abstract 

An algorithm to simulate full QCD with 3 colours at nonzero chemical 
potential on the lattice is proposed. The algorithm works for small values 
of the chemical potential and can be used to extract expectation values 
of CPT invariant operators. 



1 



1 Introduction 



Understanding the properties of matter at nonzero density and temperature 
is an essential ingredient to describe the physics of the early universe and 
the collapse of very massive stars [flj. Moreover recent heavy ion collision 
experiments offer a unique laboratory test for the predictions of the theory H . 

Effective models |], |J predict a rich structure in the temperature-chemical 
potential T—fi diagram of QCD with several phases. Numerical simulations on 
the lattice is an adequate technique to study the corresponding phase transi- 
tions. The lattice action of QCD at finite \x has been given in @, ||. By using 
Wilson fermions within a SU(N) invariant gauge theory, this action is || 



^Wilson = Pj2[ 1 -^ ReTrP ) 



flavours X 

V 

+ (A + lv ) ? x+i >C^(s)^(^)~Va 



CM = 1 + S 4u (exp (/(a/i)) - 1) , (1) 

where the first term is the pure gauge action, a is the lattice spacing, j3 is 
the inverse bare lattice coupling (/3 = 2N/go), P stands for the elementary 
plaquette, m is the fermion mass, A is the Wilson parameter and fi is the 
chemical potential, /(a/x) is an odd function of the chemical potential that 
satisfies f(x) = x + 0(x 3 ). The simplest choice is f(x) = x although others are 
possible (see for instance fTj where f(x) = arg tanhr). Notice that £, u (— A 4 ) = 
£i/(/x) _1 - Analogous expressions are valid for staggered and naive fermions. 

The expectation value of an operator O is defined by the Feynman path 
integral 

(O) = i y VU v V$ x Vil) x O exp (-Wilson) , 
Z = J VU u Vip x T>ip x exp (-Swiison) , (2) 

where terms regarding gauge fixing have been skipped as they are inessential 
for our analysis. After integrating out fermions (we assume that O can be 
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expressed in terms of gluon fields only) we have 



(O) = 1 J VU V O det D exp {-S g ) , (3) 

where S g is the pure gauge action and D is the fermion matrix 
D xy ((j) = (ma + X)5 X} y 

~\ E ( S *,v+> U t(y) ( A + 7.) W 1 + 5y, x+ „U u (x) (A - 7,) £„(//)) (4) 

Colour, flavour and spinor indices are placed where required. The theory with 
N = 3 colours and fermions in the fundamental representation in general yields 
a complex value for det D if fi ^ 0. This is a problem because importance 
sampling in the Monte Carlo integration of Eq.(|3|) cannot be applied with a 
complex weight. Several solutions have been proposed that allow to perform 
the integration in an approximate way or in particular regions of the T—y, dia- 
gram: reweighting methods |J, calculation of the canonical partition function 
from simulations at imaginary chemical potential ||, analytical continuations 
of Taylor expansions in powers of imaginary chemical potential |10| , responses 



of several observables to a nonzero small chemical potential Jl l]| and fugacity 



expansions QX2j| . Here we present another idea that should work for small values 



of the chemical potential. We will prove that if we restrict our attention to the 
calculation of expectation values of CPT invariant operators then the correct 
Boltzmann weight under the path integral is exp (—S g ) Re det D. 



2 The Boltzmann weight 



The Boltzmann weight in Eq.(£J), exp (—Sg) det D, is a complex number. How- 
ever the expectation value of any observable represented by the operator O 
in Eq.(H) must be a real quantity. Imposing that (O) be real for any observable 
O is a strong constraint on the integration measure. Clearly all possible con- 
figurations can be gathered in sets such that the contribution to the imaginary 
part of (O) from all configurations in one single set cancels out. We shall as- 
sume that all sets are formed by only two configurations and that these two are 
related by some transformation S. Firstly we want to find this transformation. 
Let C denote a thermalized arbitrary configuration on the lattice and C s the 
corresponding transformed configuration by the action of S. We require that 
(i) the value of the operator O calculated on C and C s be the same; (ii) the 
weight under the path integral for C s be the complex conjugate of the weight 
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for C and (in) C and C s have the same Haar measure, 

0[C] - 0[C s ]=0[C], 
detZ^e" 5 ^ -> detD[C s ]e-^P s ] = ( de t D[C])* e"^ [c] , 

2W„[C] -> T>U U [C S ]=VU U [C], (5) 

where we have explicitely written the dependence on the configuration. We 
will write this dependence wherever necessary. These conditions are sufficient 
to guarantee that (O) is real. 

A simple way to enforce the last constraint in (|B|) is by imposing that S is 
a discrete transformation. T>U U means Il x I\. v (iUy{x) where each single factor 
dU v (x) is the Haar measure over the gauge group. By a discrete transformation 
we mean that S does not transform each single Haar measure. In fact this would 
entail a jacobian and the search of the transformed configuration would become 
more difficult. 

The matrix D xy (n) has the property (det-D (//))* = detD(— //). This stems 
from the relation -D(/i)' = 7s-D(— /i)7sQ- Notice that this implies Re detD (Im 
det-D) is an even (odd) function of \x. Moreover it is physically clear that apply- 
ing a charge conjugation operator C changes the sign of \i. Then we expect that 
the condition det£>[C s ] = (det £>[£])* can be verified if S contains the transfor- 
mation C. However the above condition is only part of the requirements in (0). 
In particular we see from the second condition in ([5]) that the implementation 
of C must be such that the pure gluon action S g remains unaltered. 

In order to leave S g invariant, the correct transformation should involve 
some space-time rearrangement besides charge conjugation. A reasonable guess 
is the CPT transformation that we define in the following way: if our lattice 
was d dimensional and the lateral sizes were finite and equal to L 1; L 2 , 
L d (in units of lattice spacing) then the CPT transformation would act in the 
following way 

U u (x) CPT = Ut( mod(Li - X! + 1, La) + 1, 
mod(L 2 — x 2 + l, L 2 ) + 1, 

mod(L u - x„, L v ) + 1, 

j 

mod(L d - x d + 1, L d ) + 1) , (6) 

1 This is true for Wilson and naive fermions; for staggered fermions a matrix other than 
75 must be used but the result is the same. We use the euclidean definition of the gamma 
matrices such that 7„ = 7^ holds. 
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where mod(a;, L) gives the remainder on integer division of x by L. This is 
the finite volume version of U u (x) CPT = U%(—x). This transformation fulfils all 
conditions in Eq.(|5]) as a straightforward computation shows. 

As an example let us check that S g is invariant under CPT. It is enough to 
show that under the action of CPT every plaquette of the original configuration 
C maps onto one and only one plaquette of C. Let us denote P(x; is, p) the 
plaquette starting at the site x and going around through direction v and then p. 
This is P(x; u, p) = U u (x)U p (x+u)Ul(x+ p)UHx) . Applying transformation (^|) 
on each factor we obtain 



ut 


(mod(Li - 


- X\ + 1, Li) + 1, • • • 




mod(L l/ — 






mod(L p — 


x p + l,L p ) + 1,- ■ ■) x 


ul 


(mod(Li - 


- xi + 1, Li) + 1, • • • 




mod(L l/ — 






mod(L p — 


x p , L p ) + 1, • - •) x 


u u 


(mod(Li - 


- x\ + 1, Li) + 1, • • • 




mod(L^ — 






mod(L p — 


x p ,Lp) + l,---) x 


u P 


(mod(Lx - 


-Xi + l.LO + l,-.. 




mod(L I/ — 


x v + 1,L„) + 1, ■ • • 




mod(L p — 


x pi L p ) + !,-••) , 



(7) 

which, by the cyclic property of the trace, is equivalent to the plaquette 

P(mod(Li - xi + 1, Li) + 1, • • • , mod(L„ - x v , L u ) + 1, • • • , 

mod(L p - x p , Lp) + 1, • • - , mod(L d - x d + 1, L d ) + 1; u, p) (8) 

in the original configuration C. This is clearly a one-to-one mapping from 
plaquettes to plaquettes in C. This completes the proof. 

Then a configuration C with weight det D exp(— S g ) transforms under CPT 
into another configuration C CPT with weight (det D)* exp(— S g ) and the same 
Haar measure. This means that both configurations have the same odds to be 
selected by an updating algorithm. If moreover we restrict our interest only to 
CPT invariant observables O then we can say that the numerical value 0[C] ap- 
pears with a probability proportional to exp(— S^C]) (detD[C] + (det D[C})*). 
Barring a factor 2 and dispensing with the CPT transformed configuration one 
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can just count C with the weight exp(— S"g)Re det D. This is the main result of 
our paper. 

When we say CPT invariant operators we mean that the operator O is 
CPT invariant configuration by configuration. Most of the operators are CPT 
invariant after averaging over configurations (i.e. the corresponding observable 
is CPT invariant), but our constraint is stronger than this. 

Focusing our attention only on CPT invariant operators still allows to study 
several interesting problems. Indeed the average plaquette, chiral condensate, 
Polyakov loop, etc. can be calculated by evaluating expectation values of CPT 
invariant operators. On the other hand space-time valued operators 0{x) 
are in general not permitted because our method requires that the equality 
0(x) = 0{—xy holds configuration by configuration which in general is not 
true. In particular all correlation functions are not admissible. 

Are there other discrete transformations such that the imaginary part of 
the determinant is eliminated within groups of two configurations? We have 
performed the following test. We discretized a 2-dimensional gauge theoryf] on 
a 2 2 lattice and loaded it with an arbitrary configuration (an arbitrary SU(3) 
matrix on each of the 8 links). We calculated det D obtaining a complex num- 
ber. Then we rearranged the 8 SU(3) matrices in all possible ways (8!=40320) 
and for each permutation we placed them again on the 8 links and recalculated 
the determinant. Only the transformation described in @ yielded the result 
that satisfies conditions (|5|) (permutations of the 8 links that can be viewed as 
spatial or temporal translations led to the same result but they are clearly unin- 
teresting). We conclude that there are no other simple discrete transformations 
providing the complex conjugate of the original configuration weight. 



3 A note on Monte Carlo simulation 

For the importance sampling to work it is necessary that the weight be positive 
(it must behave as a probability). However Re det D can take negative values 
and this fact limits the applicability of the method. By continuity we expect 
that Re det D is positive in the majority of configurations when \i is small 
enough. This suspicion is confirmed by numerical simulation. In |TB|] results 

2 In 2 dimensions the action that introduces the chemical potential on the lattice through 
the term /Mp^ip does not present the problems discussed in Q. The density of energy obtained 
with this action on the lattice is fj 2 /2ir, the same result than in the continuum. However 
this is irrelevant for the purpose of the present test and we used the action in Eq.(pl). 
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from a simulation with 4 flavours of staggered fermions are presented. A 4 4 
lattice is used at (3 — 4.8 and am = 0.025. In Fig. 1 we show the fraction of 
configurations with a positive (negative) weight N + /N (AL/iV) as a function 
of afj, (N = N + + iV_ is the total number of configurations). It indicates that 
our method can be used at moderate values of fi. This may include the region 
in the phase diagram T-/x where present and future heavy ion experiments 
(RHIC, LHC) are going to be run (large T and small /i). 








0.1 



0.2 



0.3 



0.4 



Figure 1: Fraction of positive and negative weighted configurations as a function of 



the chemical potential. Data taken from [ 13 1 



Monte Carlo simulations with the weight Re det D would be greatly facil- 
itated if we were able to find a new matrix A such that Re det D = det A 



because then fast and well-known simulation methods for fermions |14| could 
be used. We have not found a general and efficient algorithm to construct the 
matrix A starting from D. Consequently we have to resort to algorithms which 
explicitely calculate the determinant of D. 

We shall not insist in these aspects of the problem as they will be analysed 
in a future publication containing several numerical studies p 



7 



4 Conclusions 



We have proved that the correct Boltzmann weight for updating full QCD in 
lattice simulations at finite chemical potential is 

exp(-Sg)Re det D (9) 

where S g is the pure gluon action and D is the fermion matrix. We have shown 
that this is true when we calculate expectation values of operators that are CPT 
invariant. For other operators our proof does not work. Possibly in this case 
the assumption that the imaginary part of the Boltzmann weight is eliminated 
by combining couples of configurations should be relaxed. Nonetheless the 
class of CPT invariant operators include many observables usually studied in 
the context of finite density systems. 

Our algorithm has two problems which deserves further improvement: on 
one hand the weight (|j) is not positive in general. We have shown that it is 
mostly positive for moderate values of the chemical potential. On the other 
hand the present method requires the explicit calculation of the determinant of 
the fermion matrix D which is very time consuming. In a future publication |TJ| 
we will give numerical results obtained by using Eq.(^). 
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